home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / QFLOAT / QFLTA.C < prev    next >
C/C++ Source or Header  |  1996-02-25  |  6KB  |  461 lines
  1. /*        qflta.c
  2.  * Utilities for extended precision arithmetic, called by qflt.c.
  3.  * These should all be written in machine language for speed.
  4.  *
  5.  * addm( x, y )        add significand of x to that of y
  6.  * shdn1( x )        shift significand of x down 1 bit
  7.  * shdn8( x )        shift significand of x down 8 bits
  8.  * shdn16( x )        shift significand of x down 16 bits
  9.  * shup1( x )        shift significand of x up 1 bit
  10.  * shup8( x )        shift significand of x up 8 bits
  11.  * shup16( x )        shift significand of x up 16 bits
  12.  * divm( a, b )        divide significand of a into b
  13.  * mulm( a, b )        multiply significands, result in b
  14.  * mdnorm( x )        normalize and round off
  15.  *
  16.  * Copyright (c) 1984 - 1988 by Stephen L. Moshier.  All rights reserved.
  17.  */
  18.  
  19. #include "qhead.h"
  20.  
  21. int qmovz(), cmpm(), mdnorm();
  22.  
  23. /*
  24. ;    Shift mantissa down by 1 bit
  25. */
  26.  
  27. int shdn1(x)
  28. register unsigned short *x;
  29. {
  30. register short bits;
  31. int i;
  32.  
  33. x += 2;    /* point to mantissa area */
  34.  
  35. bits = 0;
  36. for( i=0; i<NQ-1; i++ )
  37.     {
  38.     if( *x & 1 )
  39.         bits |= 1;
  40.     *x >>= 1;
  41.     if( bits & 2 )
  42.         *x |= 0x8000;
  43.     bits <<= 1;
  44.     ++x;
  45.     }    
  46. return 0;
  47. }
  48.  
  49. /*
  50. ;    Shift mantissa up by 1 bit
  51. */
  52. int shup1(x)
  53. register unsigned short *x;
  54. {
  55. register short bits;
  56. int i;
  57.  
  58. x += NQ;
  59. bits = 0;
  60.  
  61. for( i=0; i<NQ-1; i++ )
  62.     {
  63.     if( *x & 0x8000 )
  64.         bits |= 1;
  65.     *x <<= 1;
  66.     if( bits & 2 )
  67.         *x |= 1;
  68.     bits <<= 1;
  69.     --x;
  70.     }
  71. return 0;
  72. }
  73.  
  74.  
  75.  
  76. /*
  77. ;    Shift mantissa down by 8 bits
  78. */
  79.  
  80. int shdn8(x)
  81. register unsigned short *x;
  82. {
  83. register unsigned short newbyt, oldbyt;
  84. int i;
  85.  
  86. x += 2;
  87. oldbyt = 0;
  88. for( i=0; i<NQ-1; i++ )
  89.     {
  90.     newbyt = *x << 8;
  91.     *x >>= 8;
  92.     *x |= oldbyt;
  93.     oldbyt = newbyt;
  94.     ++x;
  95.     }
  96. return 0;
  97. }
  98.  
  99. /*
  100. ;    Shift mantissa up by 8 bits
  101. */
  102.  
  103. int shup8(x)
  104. register unsigned short *x;
  105. {
  106. int i;
  107. register unsigned short newbyt, oldbyt;
  108.  
  109. x += NQ;
  110. oldbyt = 0;
  111.  
  112. for( i=0; i<NQ-1; i++ )
  113.     {
  114.     newbyt = *x >> 8;
  115.     *x <<= 8;
  116.     *x |= oldbyt;
  117.     oldbyt = newbyt;
  118.     --x;
  119.     }
  120. return 0;
  121. }
  122.  
  123.  
  124.  
  125. /*
  126. ;    Shift mantissa up by 16 bits
  127. */
  128.  
  129. int shup16(x)
  130. register short *x;
  131. {
  132. int i;
  133. register short *p;
  134.  
  135. p = x+2;
  136. x += 3;
  137.  
  138. for( i=0; i<NQ-2; i++ )
  139.     *p++ = *x++;
  140.  
  141. *p = 0;
  142. return 0;
  143. }
  144.  
  145. /*
  146. ;    Shift mantissa down by 16 bits
  147. */
  148.  
  149. int shdn16(x)
  150. register unsigned short *x;
  151. {
  152. int i;
  153. register unsigned short *p;
  154.  
  155. x += NQ;
  156. p = x+1;
  157.  
  158. for( i=0; i<NQ-2; i++ )
  159.     *(--p) = *(--x);
  160.  
  161. *(--p) = 0;
  162. return 0;
  163. }
  164.  
  165.  
  166.  
  167. /*
  168. ;    add mantissas
  169. ;    x + y replaces y
  170. */
  171.  
  172. int addm( x, y )
  173. unsigned short *x, *y;
  174. {
  175. register unsigned long a;
  176. int i;
  177. unsigned int carry;
  178.  
  179. x += NQ;
  180. y += NQ;
  181. carry = 0;
  182. for( i=0; i<NQ-1; i++ )
  183.     {
  184.     a = (unsigned long )(*x) + (unsigned long )(*y) + carry;
  185.     if( a & 0x10000 )
  186.         carry = 1;
  187.     else
  188.         carry = 0;
  189.     *y = a;
  190.     --x;
  191.     --y;
  192.     }
  193. return 0;
  194. }
  195.  
  196. /*
  197. ;    subtract mantissas
  198. ;    y - x replaces y
  199. */
  200.  
  201. int subm( x, y )
  202. unsigned short *x, *y;
  203. {
  204. register unsigned long a;
  205. int i;
  206. unsigned int carry;
  207.  
  208. x += NQ;
  209. y += NQ;
  210. carry = 0;
  211. for( i=0; i<NQ-1; i++ )
  212.     {
  213.     a = (unsigned long )(*y) - (unsigned long )(*x) - carry;
  214.     if( a & 0x10000 )
  215.         carry = 1;
  216.     else
  217.         carry = 0;
  218.     *y = a;
  219.     --x;
  220.     --y;
  221.     }
  222. return 0;
  223. }
  224.  
  225.  
  226. int divm( a, ac3 )
  227. unsigned short a[], ac3[];
  228. {
  229. unsigned short act[NQ+1], ac1[NQ+1];
  230. int i, ctr, lost;
  231. unsigned short d, *p, *q, *r;
  232.  
  233. qmovz( a, ac1 );
  234. qclear( act );
  235. act[0] = ac3[0];
  236. act[1] = ac3[1];
  237. act[NQ] = 0;
  238. ac3[NQ] = 0;
  239.  
  240. /* test for word-precision denominator */
  241. for( i=4; i<NQ; i++ )
  242.     {
  243.     if( ac1[i] )
  244.         goto longdiv;
  245.     }
  246. /* Do divide with faster compare and subtract */
  247. d = ac1[3];
  248. p = &ac3[3];
  249. q = &ac3[2];
  250. r = &act[NQ];
  251. for( i=0; i<NBITS+2; i++ )
  252.     {
  253.     if( (*q != 0) || (*p >= d) )
  254.         {
  255.         *p -= d;
  256.         *q = 0;
  257.         *r = 0x8000;
  258.         }
  259.     shup1( ac3 );
  260.     shup1( act );
  261.     }
  262. goto divdon;
  263.  
  264. /* Slower compare and subtract required */
  265. longdiv:
  266. for( i=0; i<NBITS+2; i++ )
  267.     {
  268.     if( cmpm( ac3, ac1 ) >= 0 )
  269.         {
  270.         subm( ac1, ac3 );
  271.         ctr = 1;
  272.         }
  273.     else
  274.         ctr = 0;
  275.     shup1( ac3 );
  276.     shup1( act );
  277.     act[NQ-1] |= ctr;
  278.     }
  279.  
  280. divdon:
  281.  
  282. p = &ac3[2];
  283. lost = 0;
  284. for( i=2; i<NQ; i++ )
  285.     {
  286.     if( *p++ != 0 )
  287.         {
  288.         lost = 1;
  289.         break;
  290.         }
  291.     }
  292. shdn1( act );
  293. shdn1( act );
  294. act[1] -= 1;
  295. if( act[1] & 0x8000 )
  296.     act[1] = 0;
  297. mdnorm( act, lost );
  298. qmov( act, ac3 );
  299. return 0;
  300. }
  301.  
  302.  
  303.  
  304.  
  305. int mulm( b, ac3 )
  306. unsigned short b[], ac3[];
  307. {
  308. unsigned short act[NQ+1], ac2[NQ+1];
  309. unsigned short *p, *q;
  310. int ctr, nct, lost;
  311.  
  312. qmov( b, ac2 );
  313. qclear( act );
  314. act[0] = ac3[0];
  315. act[1] = ac3[1];
  316. act[NQ] = 0;
  317. /* strip trailing zero bits */
  318. nct = NBITS;
  319. p = &ac2[NQ-1];
  320. while( *p == 0 )
  321.     {
  322.     shdn16( ac2 );
  323.     nct -= 16;
  324.     }
  325. if( (*p & 0xff) == 0 )
  326.     {
  327.     shdn8( ac2 );
  328.     nct -= 8;
  329.     }
  330. lost = 0;
  331. q = &act[NQ];
  332. for( ctr=0; ctr<nct; ctr++ )
  333.     {
  334.     if( *p & 1 )
  335.         addm(ac3, act);
  336.     shdn1(ac2);
  337.     lost |= *q & 1;
  338.     shdn1(act);
  339.     }
  340. mdnorm( act, lost );
  341. qmov( act, ac3 );
  342. return 0;
  343. }
  344.  
  345.  
  346.  
  347.  
  348. int mulin( a, b )
  349. unsigned short a[], b[];
  350. {
  351. mulm(a,b);
  352. return 0;
  353. }
  354.  
  355.  
  356. unsigned short rndbit[NQ+1] = {0};
  357. short rndset = 0;
  358.  
  359. int mdnorm( x, lost )
  360. unsigned short x[];
  361. int lost;
  362. {
  363. int i;
  364. register short *p;
  365.  
  366. if( rndset == 0 )
  367.     {
  368.     qclear( rndbit );
  369.     rndbit[NQ-1] = 1;
  370.     rndbit[NQ] = 0;
  371.     rndset = 1;
  372.     }
  373. p = (short *)&x[1];
  374. for( i=0; i<3; i++ )
  375.     {
  376.     if( x[2] == 0 )
  377.         break;
  378.     shdn1(x);
  379.     *p += 1;
  380.     if( *p < 0 )
  381.         *p = 32767; 
  382.     }
  383. for( i=0; i<3; i++ )
  384.     {
  385.     if( x[3] & 0x8000 )
  386.         break;
  387.     shup1(x);
  388.     *p -= 1;
  389.     if( *p < 0 )
  390.         *p = 0;
  391.     }
  392. i = x[NQ] & 0xffff;
  393. if( i & 0x8000 )
  394.     {
  395.     if( (i == 0x8000) && (lost == 0) )
  396.         {
  397.         if( (x[NQ-1] & 1) == 0 )
  398.             goto nornd;
  399.         }
  400.     addm( rndbit, x );
  401.     }
  402. if( x[2] )
  403.     {
  404.     shdn1( x );
  405.     *p += 1;
  406.     if( *p < 0 )
  407.         *p = 32767;
  408.     }
  409. nornd:
  410. x[NQ] = 0;
  411. return 0;
  412. }
  413.  
  414. /*
  415. ;    move a to b
  416. ;
  417. ;    short a[NQ], b[NQ];
  418. ;    qmov( a, b );
  419. */
  420.  
  421. int qmov( a, b )
  422. register short *a, *b;
  423. {
  424. int i;
  425.  
  426. for( i=0; i<NQ; i++ )
  427.     *b++ = *a++;
  428. return 0;
  429. }
  430.  
  431.  
  432. int qmovz( a, b )
  433. register short *a, *b;
  434. {
  435. int i;
  436.  
  437. for( i=0; i<NQ; i++ )
  438.     *b++ = *a++;
  439. *b++ = 0;
  440. return 0;
  441. }
  442.  
  443.  
  444. /*
  445. ; Clear out entire number, including sign and exponent, pointed
  446. ; to by x
  447. ;
  448. ; short x[];
  449. ; qclear( x );
  450. */
  451.  
  452. int qclear( x )
  453. register short *x;
  454. {
  455. register int i;
  456.  
  457. for( i=0; i<NQ; i++ )
  458.     *x++ = 0;
  459. return 0;
  460. }
  461.